%
%
%--------------
%calibration.m
%--------------
%Runs within the directory that holds the calibration data
%Copy this file to the directory that contains the FIDs that would be used
%for calibrating and obtaining the trajectory
%
%The program gets a list of FID files that contains the PARSE datasets,
%checks for their lenghts (number of points), number of repetions per
%acquitions (blocks), offsets, and the gradient that was used (calib=p/calib=r).
%
%The output of this program is a the list of 'raw' real and imaginary parts,
%averaged combined FID (obtained by averaging RE and IM over number of
%repetitions), offsetlist (displacement pss0 depending upon which gradient
%is used), phi, and the k-trajectory
%-------------
%
%Rick Deshpande, Mark Bolding, Don Twieg : CDFI, UAB - August 2007
%

%{
%Get usable FIDs first. Look at the hole at the center for real vs imag
%plots. Copy paste the script below in the command window in data directory
for count=1:27
if count<10
fname=strcat('pinecone_0',num2str(count));
else
fname=strcat('pinecone_',num2str(count));
end
[re,im]=load_fid(fname);
fid=mean(re,2)+i*mean(im,2);
close all;
[np nt pss0 axis]=readprocpar(fname);
dis_str=strcat('FID: ',num2str(count),' | Offset: ',num2str(pss0),'cm | Axis: ',axis);
figure;
plot(fid);title(dis_str);
pause;
end
%}

warning off all
close all;
clear;
clc;

%% Get the directory structure
dcOffset= 20+i*20;              %Empirically computed d.c. offset
FileList=dir('pine*.*');
FileList=char(FileList.name);   % get the filenames out of struct variable
FileList=FileList(:,1:end-4);   % Take off the file type extension

% Get the number of points in each FID and the Spectral Width

[np blocks PSS0 AXIS SW]=readprocpar(FileList(3,:));

clear newlist re im PSS0 AXIS;

%% Getting the FIDs and decomposing them into Real and Imaginary components
ListSize=size(FileList,1);      % Number of SS-PARSE acquisitions in that particular study
RE = zeros(np,ListSize,blocks); % Number of Real points
IM = zeros(np,ListSize,blocks); % Number of Imaginary points

%% Load all the FIDs
OffsetList=zeros(ListSize,2);   % Number of offsets used
for numacq=1:ListSize
    [RE(:,numacq,:),IM(:,numacq,:)] =load_fid(FileList(numacq,:));
    CurrentDirectory=pwd;
    [NP NT pss0 calib]=readprocpar(FileList(numacq,:));   

    %Perform the following operations only if calib=p or calib=r, do
    %nothing if calib=n

    if calib == 'p' || calib == 'r' 

        %Save the offset list depending on the type of calib used for that
        %acquistition

        if calib == 'p'
            OffsetList(numacq,:)=[pss0 0];
        elseif calib =='r'
            OffsetList(numacq,:)=[0 pss0];
        end
        
    else
        %If calib =n, remove the FIDs and do not add anything to the offset
        %list
        
        RE(:,numacq,:)=[];
        IM(:,numacq,:)=[];
        FileList(numacq,:)=[];
        OffsetList(numacq,:)=[];
    end

    %cd(CurrentDirectory);   %Revert back to the directory where the data from current study is saved
    clear pss0 posnPSS0 posnCalib content* calib contentCalib CurrentDirectory posnGmax posnSW ;

end

%*****Run this portion of the code once you have seen the RE vs. IM plots.
%Replace the discard list with the ones whose plots don't have a well
%defined hole at the center. It discards the offsets, real and imaginary
%components of those FIDs.


UseListFID=[1 18 42];
%UseListFID=[3 4 5 11 21 26 27 28];
RE_new=zeros(size(RE,1), size(UseListFID,2), blocks);
IM_new=RE_new;
offset_new=zeros(size(UseListFID,2),2);

counter=1;
for list=1:size(RE,2)
    if ismember(list,UseListFID)==1
        RE_new(:,counter,:)=RE(:,list,:);
        IM_new(:,counter,:)=IM(:,list,:);
        offset_new(counter,:)=OffsetList(list,:);
        counter=counter+1;
    end
end

RE=RE_new;
IM=IM_new;
OffsetList=offset_new;

clear RE_new IM_new offset_new counter list

%********************************************************************

RE_avg = mean(RE,3) + real(dcOffset);   %Average RE over each displacement
IM_avg = mean(IM,3) + imag(dcOffset);   %Average IM over each displacement

FID = RE_avg + i*IM_avg;                %Combined FID = sum of real and imaginary component

phase=unwrap(angle(FID));               %Unwrapped phase information from the FIDs

clear ListSize numacq RE_avg IM_avg np

%Calibration routine

%Set the offsets (read them from the list)

xoffset=OffsetList(:,1);
yoffset=OffsetList(:,2);

%pd=diff(phase);         % find differences of phase between successive samples

%[bb,aa]=butter(8,.05);  %set up butterworth low-pass filter

for n = 1:size(OffsetList,1)
    Y(:,n) = detrend(phase(:,n));%DT070920
    Y(:,n) = Y(:,n)-Y(1,n);
end

nlength=size(Y,1);

%Set up the system and acquisition parameters

gam=4258;
twopi=2*pi;
delt=1/SW;
phi0=zeros(1,nlength-1);
% for nn=1:nlength-1
% figure(1);stem3(xoffset,yoffset,Y(nn,:));pause(0.3);
% end

for nn=1:nlength-1
    [p,S]=polyfit(xoffset,Y(nn,:)',1);
    krd0(nn)=p(1)/twopi;
    phitmp=p(2);
    %clear p S
    [p,S]=polyfit(yoffset,Y(nn,:)',1);
    kpe0(nn)=p(1)/twopi;
    phi(nn)=phitmp+p(2);
    %clear p S
end

%% show results
NN=size(krd0,2);

% clear FID S bb aa twopi xoffset yoffset phase c1 gam nn n NN pd delt Y g
% save KRDPE0907